arXiv: 1508.03133v2 [astro-ph.HE] 2 Nov 2015 


Final version of November 2, 2015 

Preprint typeset using IATgX style emulateapj v. 5/2/11 


R-PROCESS LANTHANIDE PRODUCTION AND HEATING RATES IN KILONOVAE 

Jonas Lippuner and Luke F. Roberts 1 

TAPIR, Walter Burke Institute for Theoretical Physics, California Institute of Technology, MC 350-17, 1200 E California Blvd, 

Pasadena CA 91125, USA 
Final version of November 2, 2015 

ABSTRACT 

r-Process nucleosynthesis in material ejected during neutron star mergers may lead to radioactively 
powered transients called kilonovae. The timescale and peak luminosity of these transients depend on 
the composition of the ejecta, which determines the local heating rate from nuclear decays and the 
opacity. Kasen et al. (2013, ApJ, 774, 25) and Tanaka & Hotokezaka (2013, ApJ, 775, 113) pointed 
out that lanthanides can drastically increase the opacity in these outflows. We use the new general- 
purpose nuclear reaction network SkyNet to carry out a parameter study of r-process nucleosynthesis 
for a range of initial electron fractions Ye, initial specific entropies s, and expansion timescales r. 

We find that the ejecta is lanthanide-free for Y e > 0.22 — 0.30, depending on 5 and r. The heating 
rate is insensitive to s and r, but certain, larger values of Y e lead to reduced heating rates, due to 
individual nuclides dominating the heating. We calculate approximate light curves with a simplified 
gray radiative transport scheme. The light curves peak at about a day (week) in the lanthanide-free 
(-rich) cases. The heating rate does not change much as the ejecta becomes lanthanide-free with 
increasing Y e , but the light curve peak becomes about an order of magnitude brighter because it 
peaks much earlier when the heating rate is larger. We also provide parametric fits for the heating 
rates between 0.1 and 100 days, and we provide a simple fit in Y e , s , and r to estimate whether the 
ejecta is lanthanide-rich or not. 

Subject headings: gamma-ray burst: general - gravitational waves - nuclear reactions, nucleosynthesis, 
abundances - stars: neutron 


1 . INTRODUCTION 

The merger of a compact binary system that includes 
at least one neutron star, hence the merger of a neutron 
star with a black hole (NSBH) or the merger of two neu¬ 
tron stars (NSNS), is likely to eject a significant amount 
of materia l during the final stages of coalescence (Lat- 
timer et al.||1977|) in addition to emitting gravitational 
waves that may be observed b y gravitational wave detec¬ 
tors s uch as advanced LIGO (The LIGO Scientific Col- 
laboration|2015 ) and possibly powering short gamma ray 
bursts (sGRBs) (e.g. Lee & Ramirez-Ruiz 2007 |Nakar 
20071 |Gehrels et al.|2009|). 'The material that is un E 


The material that is unbound 
during the merger "is of interest for two main reasons. 
First, the majority of the mass ejected in these events is 
very neutron-rich. Once the material decompresses from 
initial densities close to nuclear density, the large num¬ 
ber of neutrons can rapidly capture on the few heavy 
nuclides present and produce nuclei up to nuclear mass 
300. This process is called the r-process because neutrons 
are captured rapidly compared to the /3-decay timescale 
of the unstable nuclides produced by neutron capture. 
Thus the r-process quickly creates heavy, very neutron- 
rich nuclides that eventually decay back to stability after 
the neutron capture ceases (Burbidge et al. 1957). De¬ 
pending on the rate of NSBH and NSNS mergers and 
the amount of neutron-rich material ejected during these 
events, they can be the domin ant source of r-process nu¬ 
cleosynthesis in the universe (|Argast et al.||2004| [Sheri 


(A rga 

et al.|2014[ |van de Voort et al.|2015| |Ramirez-Ruiz et al. 

2015). 
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Second, observable electromagnetic signals may be as¬ 
sociated with these ejecta. A radio transient that occurs 
on a timescale of a few weeks can be powered by the inter 


action of th e ejecta with the surrounding medium (Nakar 
fc Piran| 2011 ). Additionally, radioactive decay of unsta¬ 


ble" nuclides formed during decompression of the ejecta 


can power a trai 
on a timescale < 

nsient in the optical or infrared that 

peaks 

jf a day to a week (|Li & Paczynski 

1998 

Kulkarni 2005J 

Metzger et al. 2010| |Kasen et al. 

WH 

Tanaka & Hoto 

mzaka||2013D. These are often refer: 

red to 

)rono- 

3 may 

as either “kilonovae” (Metzger et al. 2010) or “mac 

vae” (Kulkarni ! 

2005). In tact, one oFthese events 


excess in the mlrarecL after¬ 
glow of nearby GRB130603B, which was an sGRB, has 
been interpreted by some authors as a strong indicator 
of a transient powered by the decay of r-pr ocess mate¬ 


rial (Tanvir et al. 2013 Berger et al. 2013). A similar 


kilonova like excess has recently been observed in the af- 


terglow of GRB060614 ( Yang et al.|2015 Jin et al.|2015 ). 


Although almost all of the ejected material will be 
neutron-rich, there can be a significant spread in the elec¬ 
tron fraction of this neutron-rich material. The compo- 


tidally (Lattimer et al. 

1977; 

Freiburghaus et al. 

1999), 

dynamically from the region where the two neutron stars 

collide (Bauswein et al. 

2013 

Hotokezaka et al.||2013a), 


or from the accre tion disk that forms after the m erger 


(TFernandez & Metzger 2013 Perego et al. 2014 Just et al. 

2015)“ Since the material ejected by all of these mecha- 

nisms starts out as cold, catalyzed material in a neutron 
star, the final electron fraction of the material depends 
on the weak interaction timescale relative to the dynam¬ 
ical timescale of the ejecta. If the temperature and local 
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neutrino density are low, and therefore weak interactions 
are slow, the electron fraction is unaltered. This is the 
case for the tidal eject a, which is predicted t o have a very 
low electron fraction (|Korobkin et al.|[2012). Conversely, 
material ejected from the disk stays near the compact ob¬ 
ject for a long period and can achieve be ta-e quilibrium at 


lower density and higher temperature (Just et al. 2015 
Richers et al |2015 ). The dynamical ejecta from the 


con¬ 


tact region sits somewhere in between, as it is ejected 
rapidly but shocked to high temperatures and irradiated 
strongly by neutrinos, which can significantly alter the 
initia l electron fraction (Wanajo et al.|2014| Goriely et al. 
r 2015). 


~Kt low initial electron fractions (Y e < 0.2), the fi¬ 
nal composition of the ejecta is relatively insensitive to 
the initial electron fraction of the material because a 
strong r-proces s occurs and fission cycling produces a 
robust pattern ( Metzger et al.||2010 Roberts et al.||2011 


Goriely et al. 201 1 [ ). But for higher electron fractions 


(0.2 ^ Y e ^ 0.3)7"an incomplete r-process can occur and 


erties of the outflow (|Korobkin et al. 

2012 Grossman 

et al. 

2014 Kasen et al. 20l5|. In ad 

Ldition to the to- 


the ejecta at around a day—which determines the nu¬ 
clear heating rate and opacity of the material—plays a 
l arge role in determin ing the properties of the kilonova 
(Li & Paczynski 1998). Since losses due to adiabatic ex- 


pansion rob all of the initial energy from the outflow, 
almost all of the energy that powers the transient must 
come from thermalizing the products of nuclear decay 
(Metzger et al. 2010). This in turn implies that the peak 


luminosity of a kilonova is sensitive to the composition. 

The opacity of the material determines the timescale 
on which the ejecta becomes optically thin and there¬ 
fore the tim e scale on w hich the transient will peak . 


Kasen et al. (2013) and Tanaka & Hotokezaka (2013) 
have shown that continuum opacity is very sensitive to 
the presence of lanthanides, and possibly actinides, in 
the outflow. Due to their large atomic complexity, lan¬ 
thanides and actinides have a very large number of lines 
relative to iron group elements and therefore their pres¬ 
ence drastically increases the opacity of the material 
and causes predicte d kilonovae to peak on timescales 


rate does not depend as strongly on the initial electron 
fraction and it changes relatively smoothly when going 
from lanthanide-rich to lanthanide-free cases. The peak 
timescale, peak luminosity, and spectral temperature of 
our modeled kilonovae differ substantially due to the ef¬ 
fect of the lanthanides and actinides on the opacity. In 
some cases, we also find very early and bright transients 
due to a neutro n-ric h freeze-out, which wa s proposed by 
Kulkarni| (|2005|) and|Metzger et al.| (|2015j). 

In Section [2j we describe our parametrized nucleosyn¬ 
thesis calculations and discuss how lanthanide produc¬ 
tion and the nuclear heating rate varies over our cho¬ 
sen parameter space. In Section [3| we present simpli¬ 
fied kilonova light curve models and examine how these 
transients vary with outflow properties. We then con¬ 
clude in Section [U Lanthanides and actinides both have 
open /-shells and thus a similar valence electron struc- 
tnr e, w hich means their impact on the opacity is similar 
(Kasen et al. 2013). Therefore, we will use the term “lan¬ 
thanides” to refer to both lanthanides and/or actinides, 
unless otherwise noted. 

2. PARAMETERIZED EJECTA 
NUCLEOSYNTHESIS 

The details of the r-process abundance pattern, espe- 


daily the position of the third peak, 
to the nuclear mass model, reaction r 
fragment distributions that are used (e 

can be sensitive 
ates, and fission 

:.g. Goriely et al. 

2005 

Arcones & Martmez-Pinedo 2011| 

Mumpower et al. 

2012 

de Jesus Mendoza-Temis et al.| 2014; Eichler et al. 

2014 

). Here, we are less interested in the detailed final 


abundance patterns at high mass and more interested in 
the surfaces in our parameter space at which lanthanide 
production ceases. Therefore, we employ a single mass 
model and set of reaction rates. We use two models for 
fission fragments, but our main results are insensitive to 
this choice. 

Rather than po st-processing full hyd r odynamic mod¬ 
els as was done in Goriely etah I (poll) ; |Ko robkin e t al 


Wanajo et al. (|2(J14|); Just 


Martin et al.|(|2015p, we use a parametrized 
approach that allows us to systematically study the im¬ 
pact of different ejecta properties on the p roperties of 


the ejected material relevant to kilonovae. Kasen et al 


of around a week (Barnes & Kasen 2013 Tanaka & (2015) performed preliminary investigations of the elec 


Hotokezaka 2013). Older models that assumed iron- tron fraction at which lanthanide production ceases, but 


day (| 

Metzger et al. 

2010 Roberts et al. 

2011 Goriely decay heating rate and only considered a small region of 

et al. 

2011|. Significant lanthanide and acl 

tinide produc- the parameter space. 


Fernandez| (|2014|) 
the peak time of 


have suggested that measurement of 
peak time of a kilonova might provide insight into 
the composition of the outflow. 

In this work, we present a parameter study of detailed 
nucleosynthesis calculations in NSBH or NSNS merger 
scenarios and their associated kilonova light curves. We 
focus in particular on the mass fraction of lanthanides 
and actinides present in the ejecta, the radioactive heat¬ 
ing rate at 1 day, and how these properties depend on 
the initial conditions of the outflow. As expected, the 
lanthanide and actinide abundances depend strongly on 
the electron fraction, but the entropy and expansion 
timescale can also play an important role in certain 
cases. In contrast, we find that the nuclear decay heating 


the expanding material that undergoes r-process nucle¬ 
osynthesis and produces a kilonova. 

(i) The initial electron fraction Y e = N p /Nb , where 
N p is the total number of protons (free or inside nuclei) 
and Nb is the total number of baryons. We sample Y e 
uniformly between 0.01 (very neutron-rich matter) and 
0.5 (symmetric matter). We do not consider Y e > 0.5 be¬ 
cause the r-process requires a neutron-rich environment. 

(ii) The initial specific entropy s, which we sample 
logarithmically between 1 and 100 ks baryon -1 . 

(iii) The expansion timescale r, which determines 
how fast the density decreases during nuclear burning. 
We sample r logarithmically between 0.1 and 500 ms. We 
choose an analytic density profile that initially decreases 
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exponentially with time, i.e. p oc e _t / r , and then tran¬ 
sitions smoothly to a homologous, p oc t -3 , expansion. 
Requiring continuity of p and dp/dt fixes the matching 
point at t = 3r and gives 


p{t) = 


Poe */ r if t < 3r, 

"“(I) if, - 3T ' 


(i) 


where po is the initial density and e is Euler’s number. 
This parameterization of the density is chosen because 
it gives us direct control over the dynamical timescale 
at the time of r-process nucleosynthesis but still matches 
smoothly to the density profile expected for homologous 
ejecta. We have also found that this profile gives a good 
approximation to density histories of Lagrangian flui d el¬ 


ements in the ejecta of BH NS mergers simulations (Duez 


2015 Foucart et al. 


2014) 


We determine po by setting the initial temperature to 
T = 6 x 10 9 K and then finding the density for which 
nuclear statistical equilibrium (NSE) (with the given Y e ) 
produces a set of abundances that has the prescribed ini¬ 
tial entropy s. The entropy is calculated from the NSE 
distribution using a modified versio n of the Helmholtz 


equation of state (EOS) based on Timmes & Swesty 
(2000). The EOS has been modified to calculate the 
entropy for each nuclear species separately, rather than 
using average mass and charge numbers, and it also 
includes the internal partition functions of all nuclear 
species, which we obtained from the WebNucleo database 
distributee^ with REACLIB (see below). The resulting 
initial densities range from 7.1 x 10 5 to 1.4 x 10 12 gem -3 . 

Given F e , 8, and r, NSE determines po (and thus p{t)) 
and the initial abundances. We then use the newly de¬ 
veloped nuclear reaction network SkyNet for the abun¬ 
dance evolution. SkyNet is a general-purpose, modular 
nuclear reaction network that keeps track of entropy and 
temperature changes due to the nuclear reactions it is 
evolving. A detailed code description of the fu nctional¬ 


ity and featur es of SkyNet is forthcoming (Lippuner & 


Roberts|2015 , in prep.), and the source code will be pub 


licly released together with that paper. In the meantime, 
anyone who wishes to use SkyNet can contact the authors 
and request access to the code. 

We run SkyNet with nucle ar reaction rates fro m the 
JINA REACLIB databas^ 3 ] (|Cyburt et al. 2010). The 
nuclear data (masses and partition functions) were taken 
from the associated WebNucleo XML file distributed 
with REACLIB. Although REACLIB includes inverse 
rates for the strong reactions, SkyNet calculates these 
inverse rates from detailed balance, so that the rates are 
consistent with NSE. We also include different sets of 
spontaneous and neutron-induced fission rates, as REA- 


2 https://groups.nscl.msu.edu/jina/reaclib/db/library. 
php?action=viewsnapshots 

° At the time of writing, the latest REACLIB snapshot (2013-04- 
02) contains 83 incorrect /3-decay rates, which we corrected for this 
study. It appears that some lower limits of the half-lives published 
in the Nuclear Wallet Cards (http://www.nndc.bnl.gov/wallet) 
were put into REACLIB, but those lower limits can be very far 
away from realistic estimates of the half-lives. For example, REA¬ 
CLIB gives a half-life of 300 ns for 216 Pb beca use the Nuclear Wal- 
let Cards state the half-life is “> 300ns”, but|Moller et al. (2003| 
gives a half-life of about 850 s, which is much closer to the half-lives 
of similar nuclides. 


Table 1 

Parameter Values at Grid Points 






Additional 


Low-resolution points a 

high-resolution points b 

Ye 

s 

T 

Y e 

s 

T 


( k,B baryon -1 ) 

(ms) 


(kB baryon -1 ) 

(ms) 

0.01 

1.0 

0.10 

0.04 

1.3 

0.17 

0.07 

1.8 

0.29 

0.10 

2.4 

0.49 

0.13 

3.2 

0.84 

0.16 

4.2 

1.4 

0.19 

5.6 

2.4 

0.22 

7.5 

4.2 

0.25 

10 

7.1 

0.29 

13 

12 

0.32 

18 

21 

0.35 

24 

35 

0.38 

32 

59 

0.41 

42 

100 

0.44 

56 

170 

0.47 

75 

290 

0.50 

100 

500 





a The low-resolution runs of the entire parameter space use 
only these grid points. 

b For the high-resolution runs of the entire parameter space we 
double the number of grid points. The high-resolution runs 
include the grid points shown in this column in addition to the 
the same points as the low-resolution runs. 


CLIB does not presently include any fission reactions. 
There are three sets of symmetric neutron-induced fis¬ 
sion reactions: symO, sym2, and sym4, which produce 0, 
2, and 4 free neutrons, respectively, for each fission event. 
There is also a set nonsym of non-symmetric fission reac¬ 
tions that do not produce any free neutrons. Each nucle¬ 
osynthesis calculation includes one of the four neutron- 
induced fission reaction sets and the spontaneous fission 
reaction set. All the fission reactions and their ra tes ar e 
taken from the same sources used in Roberts et al. ( 2011). 
We use be t a-deca y and electron cap ture rates from 


Fuller et al. 


(1982), Oda et al. (1994) and Langanke 


& lVlartmez-Hmedo| (2000) whenever they are available 
For nuclei for which these rates are not available, the 
effects of electron blocking and positron capture are ap¬ 
proximately included by assuming that only a ground 
state to g round stat e transition occurs as described in 


Arcones et al. (2010). These rates are then normalized 
such that they are equal to the vacuum decay rates given 
in REACLIB at low temperature and density, which 
can be thought of as setting the effective matrix ele¬ 
ment for the ground state to ground state transition. 
Because this procedure assumes a maximal Q-value for 
these weak rates, this provides a lower limit on the effect 
of the surrounding medium on the combined beta-decay 
and lepton capture rate. For this study, we run SkyNet 
with 7843 nuclear species, ranging up to Z = 112 and 
A = 337, and 110,793 nuclear reactions. 


2.1. Parameter space 

We use a 9 x 9 x 9 grid to cover the entire parameter 
space and run SkyNet for each point with all four sets 
of neutron-induced fission reactions (symO, sym2, sym4, 
nonsym). We also run the symO fission reactions with 
a finer 17x17x17 grid. The parameter values at the 
grid points are shown in Table [1] The different fission 
reactions only result in small quantitative and no quali- 
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Figure 1. The final abundances of some selected nucleosynthesis calculations. Left: Y e = 0.01,0.19,0.25,0.50, s = 10 ks baryon -1 , and 
r = 7.1 ms. The full r-process is made, with substantial amounts of lanthanides and actinides, for Y e = 0.01 and Y e = 0.19. The Y e = 0.25 
trajectory is neutron-rich enough to make the second r-process peak, but not the third and not a significant amount of lanthanides. In 
the symmetric case ( Y e = 0.5), mostly 4 He and iron-peak elements are produced. Right: Y e = 0.25, s = 1.0, 3.2,10,100 ks baryon -1 , and 
t = 7.1 ms. With s = 1 ks baryon -1 a jagged r-process is obtained because there are only few free neutrons per seed nucleus available and 
nuclides with even neutron numbers are favored. Even though there are not many free neutrons available, there is still a significant amount 
of lanthanides in the s = Iks baryon -1 case because the initial seed nuclei are very heavy. At higher entropies, the initial seeds become 
lighter and the initial free neutron abundance increases. However, the increase in the initial free neutron abundance is not enough to offset 
the decrease in the initial mass of the seeds and so we obtain a less complete r-process. The situation is reversed at s = 100 ks baryon -1 , 
where there is a very high neutron-to-seed ratio. In that case, a significant fraction of a particles are also captured on the seed nuclei. This 
leads to a full r-process in the s = 100 ks baryon -1 case. 


Temperature = 3.65E+08 K 



Unstable nuclide 
□ Stable nuclide 
■ Missing nuclide 
ii Closed neutron shell 
= Closed proton shell 
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Figure 2. A frame from the animation of the nucleosynthesis calculation for Y e = 0.01, s = 10 ks baryon -1 , and r = 7.1ms. The frame 
shows the full extent of the r-process just when free neutrons get exhausted. The plot in the upper left corner shows the temperature, 
density, and heating rate as function of time. The colored bands in the chart of nuclides correspond to the mass bins in the histogram at 
the bottom. The histogram shows the mass fractions on a linear scale while the blue curve shows the abundances as a function of mass on 
a logarithmic scale. The full animations are available at http://stellarcollapse.org/lippunerroberts2015 
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tative differences. Thus we only discuss and show plots 
of the high-resolution symO runs. Finally, we carry out a 
set of runs with high Y e resolution (A Y e = 0.005 result¬ 
ing in 99 Y e points) for s = 1,10, 30,100 ks baryon -1 
and r = 0.1,1,10 ms with the symO fission reactions. 
The data underlying all the results shown and discussed 
here (nucleosynthesis results, heating rate fit coefficients, 
light curve model results, and integrated fractional heat¬ 
ing contributions of all nuclides) are available at http: 
//stellarcollapse.org/lippunerroberts2015. 

Figure [l] shows the final abundances of a few selected 
cases, which span the whole range of Y e and s at inter¬ 
mediate values of the other two parameters. For the 
s = 10 k,B baryon -1 and r = 7.1ms trajectories (left 
panel of Figure [Tl ), the full r-process up to the third 
peak (A ~ 190) tor Y e = 0.01 and Y e = 0.19 is pro¬ 
duced. We note good agreement of the second, third, and 
rare-earth peak positions with the solar r-process abun¬ 
dances, although the third peak is slightly overproduced 
relative to the second peak. The abundance patterns of 
Y e = 0.01 and Y e = 0.19 are very similar because both 
cases are neutron-rich enough to produce nuclides with 
A > 250, which eventually undergo fission. As the ejecta 
becomes less neutron-rich (Y e = 0.25 and Y e = 0.50), 
the full r-process is no longer produced; there are not 
enough neutrons available per seed nucleus to reach the 
third peak. At Y e = 0.25, the first and second r-process 
peaks are produced. The right panel of Figure [l] shows 
the final abundances of cases with Y e = 0.25, r = 7.1 ms, 
and different initial entropies. Here, the electron fraction 
is too high to get to the third r-process peak at most 
entropies (all the cases with entropies between 10 and 
75 ks baryon -1 have virtually identical final abundances 
as the s = 10 ks baryon -1 case). At s = 100 ks baryon -1 
the third r-process peak is obtained because the initial 
composition contains few seed nuclei and alpha particles 
are unable to efficiently combine to produce seed nuclei. 
Thus, the neutron-to-seed ratio is significantly enhanced. 

Animations of the full nucleosynthesis calculations for 
all seven cases shown in Figure [l] are available at http: 
//stellarcollapse.org/lippunerroberts2015. Fig¬ 
ure [2] shows a frame from one of the animations. 

2.2. Lanthanide turnoff and heating rate as a 
function ofY e 

Figure [3] shows the final lanthanide and actinide mass 
fractions Xpa and X Ac , respectively, along with the neu¬ 
tron mass fraction X n at 10 minutes, which is the mean 
lifetime of a free neutron. Also shown is Afi n , which is 
an estimate of the final average mass number A of the 
material. It is defined as 


^fin — 


•^seed(O) T Fa(0)/18 


( 2 ) 


where Y^(0) is the initial a-particle abundance and 
^seed(O) is the initial seed abundance (sum of abundances 
of all nuclides with A > 12). Since the ce-process ceases 
around K r in neutron rich conditions ( |Woosley fc Hoff- 


man 


1992), it takes around eighteen a particles to make a 


seed nucleus. Therefore, the quantity in the denominator 
of Equation pi) is approximately the number abundance 
of heavy nuclei present at the end of the r-process. We 
then arrive at Equation © by assuming that the total 


Table 2 

Afi n and Y e at Lanthanide and Actinide Turnoff 


Lanthanide turnoff a Actinide turnoff 3. 


s 

( k,B baryon - 

- 1 ) (ms) 

Ye 

Afi n 

Ye 

Afin 

1.0 

0.1 

0.27 

94 

0.25 

123 

1.0 

1 

0.28 

91 

0.24 

137 

1.0 

10 

0.28 

93 

0.18 

192 

1.8 

0.1 

0.25 

106 

0.21 

123 

1.8 

1 

0.27 

100 

0.21 

125 

1.8 

10 

0.27 

98 

0.17 

170 

3.0 

0.1 

0.23 

118 

0.20 

135 

3.0 

1 

0.25 

111 

0.21 

130 

3.0 

10 

0.27 

106 

0.18 

150 

5.6 

0.1 

0.22 

135 

0.14 

196 

5.6 

1 

0.23 

127 

0.21 

138 

5.6 

10 

0.24 

124 

0.21 

140 

10 

0.1 

0.13 

223 


- 

10 

1 

0.24 

121 

0.21 

139 

10 

10 

0.24 

120 

0.21 

139 

18 

0.1 





18 

1 

0.2 1 

102 

0.20 

12.0 

18 

10 

0.24 

102 

0.21 

125 

30 

0.1 


— 



30 

1 

0.2 1 

93 

0.18 

12,2 

30 

10 

0.24 

93 

0.20 

113 

56 

0.1 





56 

1 

0.2 1 

0 1 

0.10) 

1 12, 

56 

10 

0.24 

94 

0.21 

109 

100 

0.1 





100 

1 

0.28 

0 1 

0.18 

1 18 

100 

10 

0.29 

92 

0.26 

102 


a Turnoff is when the mass fraction Ap a or X Ac drops below 
lO -3 . The columns show the maximum Y e and corresponding 
minimum for which Xi > 10 —3 . A dash (—) denotes that 
Xi < 10 -3 for all Ye, which means there is a neutron-rich 
freeze-out. 


mass fraction of heavy nuclei at the end of the calculation 
is unity. Clearly, this assumption breaks down if there is 
fission cycling, because then the number of seeds at the 
end is much larger than the number of initial seeds plus 
those produced by the ce-process. However, we are inter¬ 
ested in the value of Afi n at the actinide and lanthanide 
turnoff, which preclude significant fission cycling because 
fission cycling only happens if nuclides heavier than ac¬ 
tinides are produced, and so there is no problem in using 
the definition in Equation ©• At low electron fractions, 
<a-rich freeze-out does not occur due to the low initial 
abundance of a particles. We emphasize that Afi n only 
depends on the initial abundances, and thus it is useful 
to determine whether a certain trajectory is likely to pro¬ 
duce large quantities of lanthanides or actinides, without 
having to perform any nucleosynthesis calculation. 

Table [2] shows the values of Y e and Afi n at which 
lanthanide and actinide production ceases (mass frac¬ 
tion goes_ below 10 -3 ). In other words, if Y e is lower 
than or A& n larger than what is shown in Table [2j then 
X\, a _> 10 -3 or X Ac > 10 -3 . The lanthanide turnoff is 
at Afi n ^ 100 and the actinide turnoff is at A& n ~ 130. 
The cases where X La < 10 -3 or A Ac < 10 -3 for all Y e 
are denoted by ” in Table |2j and they correspond to 
the strong neutron-rich freeze-outs in Figure [3j which 
means that the r-process did not happen (or at least 
not efficiently) in those cases because after about 10 min 
we are just left with free neutrons that will now de¬ 
cay to protons. In the case slOrO.l (which stands for 
s = 10 %b baryon -1 and r = 0.1ms) where lanthanides 
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Figure 3. Results of the high-resolution Y e runs. The lanthanide and actinide mass fractions, X L a and Xa Cj and their sum, X^a+Acj 
are fairly constant up to some critical value of Y e in most cases because of fission cycling. The neutron abundance X n at 10 minutes (the 
mean lifetime of a free neutron) is an indicator for a neutron-rich freeze-out, which occurs at high initial entropies and short expansion 
timescales, where the neutrons do not have time to capture on the seed nuclei. The heating rate Me at 1 day with M = 10 —2 M@ is fairly 
insensitive to Y e: except at Ifigh electron fractions ( Y e > 0.4) where some individual nuclides start to dominate the heating. The estimated 
final average mass number falls off monotonically with Y e in all cases except s = 100 ks baryon -1 , where it rebounds at Y e very 
close to 0.5. There, the number of seed nuclei decreases drastically because a-particles are initially produced in higher quantities, which 
increases the neutron-to-seed ratio. In those cases, the predicted number of fission cycles Nf is artificially increased at high Y e , because of 
production of seed nuclei by the triple-a process. Where equation [3] accurately predicts the number of fission cycles, Nf falls off rapidly 
with Y e and the point where it becomes zero is correlated with the actinide turnoff, because actinides are at the low end of the fissionable 
material mass range. Note that we plot and Nf on_ linear scales rather than log scales as all the other quantities. Also, we added a 
negative offset of 5 to both and Nf and we scaled A by 1/100 so that they fit onto our left vertical axis. 
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are made, but no actinides above a mass fraction of 10 -3 , 
we see a weaker neutron-rich freeze-out in Figure [3j The 
neutron-rich freeze-outs happen at high initial entropies 
and short expansion timescales, where the ejecta is very 
hot and expands quickly, which leaves little time for 
neutrons to capture on seed nuclides. There is also a 
neutron-rich freeze-out in s30rl and slOOrl models, but 
the freeze-out is weak enough to allow lanthanides and 
actinides to be pro duced, albeit in lower quantities. Met¬ 
zger et al. (2015) suggested that a kilonova containing 
some mass with such short dynamical timescales could 
be preceded by an ultraviolet transient powered by these 
frozen-out neutrons. 

Figure [3] shows that the heating rate from decay at 
1 day is quite insensitive to Y e at Y e < 0.35 and also 
fairly insensitive to the amount of lanthanides and ac¬ 
tinides produced. As long as Xl & +Ac is more or less 
constant as a function of Y e , Me at 1 day is also fairly 
constant. When the lanthanides turn off, there is a small 
bump in the heating rate in most cases and at larger Y e , 
after lanthanides have completely gone away, the heat¬ 
ing rate drops only slightly (an order of magnitude or 
less). One might expect a larger decline of the heating 
rate once the full r-process stops happening, because the 
material is less neutron-rich overall, more stable nuclei 
are produced directly, and thus the total radioactive de¬ 
cay energy should be lower. This is indeed true and we 
verified it by looking at the integrated nuclear heating 
amount as a function of Y e (for fixed 5 and r). We find 
that in most cases the total amount of heating drops by 
1.5 to 2 orders of magnitude as Y e goes from low values to 
high values. There is a smaller drop in the heating rates 
shown in Figure [3j because there we only plot the instan¬ 
taneous heating rate at 1 day. Since the /3-decay energy is 
correlated with the decay timescale, we always see a sim¬ 
ilar instantaneous decay rate at the same point in time, 
as long as we have a collection of nuclides with half-lives 
at around a day. The picture changes at Y e > 0.35 be¬ 
cause there the final composition is dominated by one or 
a few individual nuclides, as opposed to a large ensemble 
of nuclides, which then determi ne th e heating rate. This 
is discussed in detail in Section [2~4l 

Since our parameter space is three-dimensional, we can 
go beyond giving a simple Y e cutoff for lanthanide pro¬ 
duction. We use a heuristic method to fit for the co¬ 
efficients of three inequalities in Y e , Ins, and lnr that 
separate the lanthanide-rich and lanthanide-free regions 
of the parameter space. We find that 

Wa+Ac > 10 -3 if and only if 
-1.00 Y e - 0.00744 In s kB + 0.000638 In r ms + 0.259 > 0 
and 


-0.990 Y e -h 0.117 In s kB - 0.0783 In r ms + 0.452 > 0 

and 

—0.799 Y e — 0.288 In Sk B + 0.528 In r ms + 1.88 > 0, 

where Sk B is the entropy s in units of ks baryon -1 and 
r ms is the expansion timescale r in units of milliseconds. 
The above statement only fails for 97 out of 4913 points 
in our parameter space, i.e. it is true for 98% of the 
parameter space. Most of the points where the above 
fails are very close to one of the planes, but there are 
a few points further away from the boundaries that fail 
too. Those points are all at very low F e , high entropy, 


and very short expansion timescale, where we get strong 
neutron-rich freeze-out. The results of t he f ull parameter 
space are discussed in detail in Section [2~4l 


2.3. Fission cycling 

If the r-process is strong enough to produce nuclides 
with masses near 300, these nuclides fission and the fis¬ 
sion products then capture more neutrons, eventually 
getting up to A ~ 300 and fissioning again, creating a 
fission cycle. Thus fission cycling limits the maximum 
mass of nuclides produced in the r-process, which washes 
out the initial conditions of the ejecta and hence the fi¬ 
nal abundances are determined by nuclear physics rather 
than the properties of the outflow. 

The quantity Nf shown in Figure [3] is an estimate for 
the number of fission cycles that occurred during nucle¬ 
osynthesis. It is defined as 


AT _ ^seed(^ — ^n) 

f ~ Y seed (* = o) “ ’ 


(3) 


where F seed (£ m t n ) is the abundance of all seed nuclides 
(A > 12) at the time that neutrons are exhausted (when 
X n < 10 -4 ) and Y seed (t = 0) is the initial abundance of 
seed nuclei. This estimate for the number of fission cy¬ 
cles rests on the assumption that only fission can create 
additional seed nuclides. When a neutron captures on a 
seed nuclide, it creates a heavier nuclide, but it will not 
increase the total number (and hence abundance) of seed 
nuclides in the ejecta. However, if a heavy nuclide (which 
is counted as a seed nuclide) fissions, then there are two 
seed nuclides in its place. Thus comparing the number of 
heavy nuclides at the time when neutron capture ceases 
to the initial number of heavy nuclides tells us how many 
additional heavy nuclides were produced. For example, 
if F seed (£ = t n ) m Y seed (t = 0), then no additional heavy 
nuclides were produced and thus there was no fission cy¬ 
cling, hence N f = 0. But if F seed (t = t n ) = 3Y seed (t = 0), 
for example, then (on average) each initial heavy nuclide 
produced two additional heavy nuclides and so there were 
two fission cycles, hence Nf = 2. Note that this method 
of estimating the number of fission cycles breaks down 
if nuclides with A > 12 are produced from nuclides with 
A < 12, e.g. 12 C from three 4 He. This happens most 
prominently at Y e close to 0.5 and at high entropies, 
where fission will clearly not occur. 

As expected, there are many fission cycles at low Y e 
where large amounts of lanthanides and actinides are 
produced. In the regions with significant fission cycling, 
Aha 5 Xac, and e are fairly insensitive to Y e because fission 
cycling effectively limits the maximum mass of nuclides 
that are produced to A ~ 300. As the ejecta becomes less 
neutron-rich, fewer fission cycles occur because there are 
not enough free neutrons to produce fissionable material 
with A > 250. 

In most panels in Figure [3] we see that the production 
of actinides is closely tied to fission cycling; actinides 
go away just after fission cycling stops. If the r-process 
cannot get to 4 ^ 250, it cannot create actinides and 
it cannot create fissionable material. Furthermore, in 
most panels, but especially in slrl and sir 10 there is 
an increase in Aa c and decrease in Al & at the electron 
fraction where fission cycling stops and just before ac¬ 
tinides are not produced. Just as fission cycling stops, 
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Figure 4. Slices of constant electron fraction showing the lanthanide and actinide mass fraction XL a+ Ac and the heating rate Me at 
1 day with M = 10 -2 M©. For Y e = 0.01, the high-s/small-r corner is lanthanide-free because the high entropy produces very light seed 
nuclides, fewer seed nuclei are produced due to an a-rich freeze-out, and neutron capture begins at low density due to the high entropy 
(see the text for more discussion). The low-s/large-r corner is lanthanide-free because the slow expansion timescale results in significant 
late-time heating, which drives the ejecta back to NSE, but at those late times, /3-decays have significantly raised the electron fraction and 
so the r-process starts again but at a much higher Y e , which does not produce lanthanides. The Y e = 0.25 slice is the transition between 
lanthanide-rich and lanthanide-free. At low entropies we can still make significant amounts of lanthanides because the seed nuclides are 
heavy, and at very high entropies we initially have a lot of free neutrons and a particles, which can produce significant amounts of heavy 
elements. Finally, at Y e = 0.50 the material is simply not neutron-rich enough to make any lanthanides. The heating rate at 1 day is 
quite insensitive to s and r, except at low Y e , where it is significantly smaller at high entropies and fast expansion timescales because a 
neutron-rich freeze-out happens. The uniformity in the heating rate is due to the fact that there is an ensemble of nuclides contributing to 
the heating. And since we are considering the heating at 1 day, we tend to pick up nuclides with similar decay energies (because the decay 
energy is correlated with the half-life), leading to similar heating rates even if the composition varies. 


the r-process can get to about A = 250, but not much 
above. This means that actinides can still be produced, 
but they are not being fissioned (because only lighter ac¬ 
tinides are produced or there are no more free neutrons 
to initiate fission). Lanthanides have a mass around 150 
and so they can be created from fission products. When 
fission is just turning off, we lose a small source of lan¬ 
thanides leading to the (small) decline in X L a that can 
be prominently seen in sir 10 in Figure [3] at Y e = 0.17. 

2.4. Lanthanide production and heating rate in the full 
parameter space 

Since the amount of lanthanides determines the opac¬ 
ity of the ejecta and the nuclear heating rate determines 
the amount of energy available for the electromagnetic 
transient, we are especially interested in how these two 
quantities are correlated in our parameter space. Fig¬ 
ures atom show slices of the final lanthanide and ac¬ 
tinide mass fractions, Xl^+Acj and heating rates at 1 day 
for the extreme and intermediate values of Y e , s, and 


r. All the other slice plots are available at http:// 
stellarcollapse.org/lippunerroberts2015. In the 
following, the term “lanthanide” will stand for both 
lanthanides and actinides, unless actinides are specifi¬ 
cally mentioned. Unsurprisingly, AAa+Ac depends most 
strongly on Y e and the ejecta is lanthanide-free for Y e > 
0.26. However, even for a very low Y e of 0.01, there are 
some combinations of s and r that yield a lanthanide-free 
ejecta (see upper left panel of Figure [4]). Specifically, at 
high entropies (s > 20 baryon -1 ) and small expan¬ 
sion timescales (r < 1ms), no lanthanides are produced. 
The reason for this is that neutron capture begins at 
a lower density because of the high entropy (for a fixed 
temperature at which neutron capture begins) and there¬ 
fore the neutron capture timescale is increased. This—in 
combination with light seed nuclei, a large initial neu¬ 
tron abundance, a potentially <a-rich freeze-out, and a 
short dynamical timescale—prevents production of lan¬ 
thanides and sometimes results in a neutron-rich freeze- 
out. At lower entropies, the seed nuclei are heavier and 



























Lanthanides and Heating Rates in Kilonovae 9 



Figure 5. Slices of constant entropy showing the lanthanide and actinide mass fraction AT a+ Ac and the heating rate Me at 1 day with 
M = 10 -2 Mq. At s = 1 k,B baryon -1 , no lanthanides are produced at large expansion timescales because the material heats up significantly 
at late times, which restarts the r-process at late times after Y e has risen to about 0.3. At s = 100 ks baryon -1 , no lanthanides are produced 
when the dynamical timescale is short for the reasons discussed in the caption of figure^] In all cases, there is a critical value of Y e where 
lanthanide production abruptly ceases. The heating rate at 1 day only shows some structure at high Y e where certain individual nuclides 
dominate the heating. The reduced heating in the low-Ye/small-r corner of s = 100 ks baryon -1 is due to a neutron-rich freeze-out that 
occurs there. 


the density is higher during the neutron capture period, 
allowing neutrons to capture on them even at small ex¬ 
pansion timescales. And at larger expansion timescales, 
there is more time for the neutrons to capture on the 
light seed nuclei even at very high entropies. This is 
reflected in the upper right panel of Figure [5] where no 
lanthanides are produced at small expansion timescales 
at s = 100 ks baryon -1 , and in the upper left panel of 
Figure [6] where no lanthanides are produced at high en¬ 
tropies at r = 0.1 ms. 

There is another lanthanide-free corner in the upper 
left panel of Figure 0] at very large expansion timescales 
(r > 400 ms) and low entropies (s < 3 ks baryon -1 ). 
Here, the full r-process is being made, since the ma¬ 
terial is very neutron-rich, but because the expansion 
timescale is so long, the density is still quite high (about 
10 10 gem -3 ) when neutron burning ceases. All the heavy 
elements then decay and considerably heat up the mate¬ 
rial (to above 7GK), which destroys all heavy nuclides 
via photodissociation and brings the composition back to 
NSE. Only after tens of seconds has the material cooled 
down enough for neutron captures to happen again, but 
by then, /3-decays have raised Y e to about 0.3. Thus we 
now get an r-process with an initial Y e of 0.3, which is not 
neutron-rich enough to produce lanthanides. At faster 


expansion rates (smaller r) the density falls off faster, 
resulting in less dramatic heating that cannot force the 
composition into NSE. Because we obtain the initial den¬ 
sity from solving for NSE at the prescribed entropy, Y e , 
and T = 6 GK, the initial density is lower at higher en¬ 
tropies (s > 3 ks baryon -1 ) and so even though the den¬ 
sity remains close to the initial value at r = 500 ms, the 
density is not high enough to produce heating that re¬ 
sults in NSE. This is reflected in the upper left panel 
of Figure [5] where the ejecta is lanthanide-free at large 
expansion timescales at s = 1 ks baryon -1 , and in the 
upper right panel of Figure [6] where no lanthanides are 
produced at low entropies at r = 500 ms. 

The Y e = 0.25 slice in Figure [4] is right at the transi¬ 
tion from lanthanide-rich to lanthanide-free ejecta. The 
upper panels of Figures [5] and [6] show clearly that this 
transition is very sharp at K ^ 0.22 — 0.30. In the up¬ 
per middle panel of Figure |4j the low-s/large-r corner 
that is lanthanide-free has expanded and so has the high¬ 
s/small-r corner, relative to the Y e = 0.01 panel. Addi¬ 
tionally, lanthanide production is suppressed at interme¬ 
diate entropies (5 /^baryon -1 < s < 90 Jzb baryon -1 ). 
At low entropies, we still get an r-process because the 
seed nuclei are very heavy and thus require fewer neu¬ 
trons to capture on them to make the r-process distri- 
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Figure 6. Slices of constant expansion timescale showing the lanthanide and actinide mass fraction XL a+ Ac an< l the heating rate Me at 
1 day with M = 10 -2 Mq. At r = 0.10 ms, there are no lanthanides at high entropies because the neutrons have no time to capture on the 
light seed nuclides. At r = 500 ms, there are no lanthanides at low entropies because the heavy, neutron-rich seed nuclei lead to substantial 
late-time heating, which restarts the r-process at Y e ~ 0.3, which is not neutron-rich enough to produce lanthanides. In all cases, there is a 
fairly uniform lanthanide cutoff as Y e goes beyond a critical value. The heating rate at 1 day only shows structure at high Y e where certain 
individual nuclides dominate the heating. 


bution. At very high entropies, the initial composition 
includes a large fraction of free neutrons and a particles. 
At high entropies, production of seed nuclei via neutron 
catalyzed triple-ce is suppressed (Hoffman et al. 1997), 
which reduces the number of seed nuclei and thereby in- 
creases the neutron-to-seed ratio. These conditions allow 
for the production of the r-process nuclei. With Y e > 0.3, 
lanthanides are not produced at any entropy and ex¬ 
pansion timescale combination, sinc e the ejecta is not 
neutron-rich enough. In Section |2.2| we discussed in de¬ 
tail how the final lanthanide and actinide mass fractions 
depend on Y e . 

The lower row of panels in Figures [4] to [6] shows the 
heating rate (actually Me where M = 10 -2 M©) at 1 day. 
For 0.04 < Y e < 0.35 all the Y e slices are very similar to 
the lower middle panel of Figure [4j with virtually no 
structure. At Y e = 0.01, the high-s/small-r corner has 
significantly less heating because the initial density is 
very low (po ^ 8 x 10 5 g cm -3 ) and this, coupled with the 
rapid expansion timescale (r = 0.1 ms) and the fact that 
the initial composition contains few seed nuclei (98% of 
the mass is neutrons), means there is little opportunity 
for neutron capture. For larger expansion timescales, the 
initial conditions remain the same (low initial density and 
98% of the mass is neutrons), but because the density de¬ 
creases more slowly, there is sufficient time for neutrons 


to capture on the few seed nuclei available and make a full 
r-process. At lower initial entropies, the initial density 
is larger (e.g. 4 x 10 6 gcm -3 at s = 32 baryon -1 ) so 
that the density remains higher even with a rapid expan¬ 
sion, giving the neutrons a better chance to capture on 
seed nuclei—of which there are slightly more available— 
leading to a moderate r-process. This is reflected in the 
low-F e /small-r corner of the lower right panel in Figure [5] 
and in the low-F e /high-s corner of the lower left panel in 
Figure [6] 

For Y e > 0.35 we start to see large variations in the 
heating rate at 1 day as a function of T e , which can be 
seen in all lower panels in Figures [5] and [6] But the heat¬ 
ing is still quite insensitive to s and r, as the lower right 
panel of Figure [4] shows. This variation as a function of 
Y e at high Y e can also be seen in Figure |3j There is a pro¬ 
nounced peak in the heating rate at 1 day at Y e = 0.425 
in all but the s = 100 ks baryon -1 cases. This peak is 
due to the decay of 66 Cu (half-life of 5 minutes) which 
comes from the decay of 66 Ni, which has a half-life of 
55 hours. 66 Ni has 28 protons and 38 neutrons and so 
its electron fraction is 28/66 « 0.424, which is very close 
to Y e = 0.425, the initial electron fraction of the mate¬ 
rial. Thus the initial NSE distribution contains a larger 
quantity of 66 Ni at Y e = 0.425 than at different Y e , which 
leads to excessive heating via the decay chain described 
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above because 66 Cu has a fairly large Q-value of 2.6 MeV. 
At s = 100 ks baryon -1 the initial neutron-to-seed ratio 
is much larger than at lower entropies and so significant 
neutron burning occurs even at high F e , which washes 
out the strong dependence of the heating rate at 1 day 
on Y e . 

In Figure [3j there are also large minima in the heating 
rate at 1 day in all but the s = 100 ks baryon -1 cases 
at electron fractions between 0.45 and 0.48, depending 
on s and r. These minima can also be seen in Fig¬ 
ures [5] and [6j In those cases, NSE preferentially pro¬ 
duces stable isotopes in the initial composition, which 
drastically reduces the heating. For example, the cases 
with s = l&gbaryon -1 have the minima at Y e = 0.465 
and over 80% of the initial mass is either stable or has a 
half-life of more than 100 days. The most abundant nu¬ 
clide (37% of the mass) is 56 Fe, which is stable and has 
Y e = 26/56 ~ 0.464, which is why the minimum occurs 
at Y e = 0.465, because that favors 56 Fe the most. As 
another example, the s = 10 ks baryon -1 cases have the 
minima at Y e = 0.45, where 58 Fe and 62 Ni are preferen¬ 
tially produced by NSE, which have electron fractions of 
0.448 and 0.452, respectively. 

As in Section |2.2[ we do not find a significant corre¬ 
lation between the amount of lanthanides and actinides 
produced with the heating rate at 1 day. The heating 
rate at 1 day is very uniform at values of Y e where lan¬ 
thanides are produced. Since we are looking at the heat¬ 
ing rate at a specific time, we will always pick out the 
nuclides with a half-life of about 1 day (or decay products 
of nuclides that decay on a one-day timescale). Because 
the decay energy is correlated with the half-life and be¬ 
cause we always have a collection of different nuclides, 
we obtain roughly the same heating rate at 1 day regard¬ 
less of the exact composition of the ejecta. This is no 
longer true at higher Y ej where the composition can be 
dominated by individual nuclides, which then determine 
the heating rate. 

2.5. Fitted nuclear heating rates 

For each nucleosynthesis calculation, we calculate a 
parametric fit for the nuclear heating rate e(t) between 
0.1 and 100days (the fit window). The fit has the form 

e(t) = At~ a + B ie ~ t/fil + B 2 e~ t/ ^ + B 3 e~ t/03 , (4) 

where t and /?, are in days and e(t) is in ergs -1 g -1 . We 
use at most six parameters for the fit, so either A and 
a are zero or one or more of and Pi are zero. We 
use a weighted fit where the range 0.1 to 100 days has a 
weight of one and the weight decreases linearly to zero in 
logspace from 0.1 to 0.05 days and from 100 to 200 days. 
We use a heuristic method to find the global best fit for 
all six types of fits (power law with 0, 1, or 2 exponentials, 
or 1, 2, or 3 exponentials without a power law term). The 
best of these six fits is then selected with a small penalty 
term for the number of parameter pairs. The fitting error 
is multiplied by 1.1 for each parameter pair in excess of 
one, so that we do not pick up meaningless parameters 
that improve the fit by less than 10%. 

For consistency, we calculate the fitting error at the 
same times t{ for all cases and we interpolate the actual 
heating rate to those times, which are 500 points uni¬ 
formly sampled in logspace between 10“ 2 and 10 3 days 



Time [day] 


Figure 7. Some heating rate fits showing the fits with the largest 
and smallest error, and fits with errors in between. The heating 
rate is only fitted inside the fit window (0.1 to 100days). We use a 
power law with up to two exponential terms, or up to three expo¬ 
nential terms without a power law show in Equation |4j) , whichever 
produces the best fit. The fit error (A In e/ In e) is denned in Equa¬ 
tion § As the second and third case from the top show, the fit 
can be quite bad outside the fit window. This is no surprise since 
we do not fit the data outside the fit window and because we only 
use up to three exponential terms. In reality, there are hundreds of 
individual nuclides contributing to the total heating rate and each 
one contributes a different exponential term. 


(however, points before 0.05 days and after 200 days have 
zero weight and thus do not contribute to the fitting er¬ 
ror, as explained above). The fit error used for finding 
the optimal fit parameters is the sum of squares of the 
log difference, i.e. 

fit error = Wj (lne(^) — lne(t^)) 2 , (5) 

i 

where Wi is the weight of time ti. This error measure 
works well for the optimization algorithm to find the best 
parameters, but it carries little physical meaning. To be 
able to intuitively judge the quality of a particular fit, 
we define the mean fractional log error as 

/ A In e \ / | In e(£») — In e(t») | \ 

\i^r/ = \——/• (6) 

where the average runs over all times ti such that 
0.1 days < U < 100 days. We only fit the total heating 
rate, but we also provide the average heating contribu¬ 
tion due to fission reactions in the fit window. 

The best and worst heating rate fits, as well as some 
fits of intermediate quality, are shown in Figure [7] About 
80% of all high-resolution symO fits have (A In e/ In e> < 
0.5% and about 95% have a mean fractional log error of 
at most 1%. Since we do not include /^-delayed fission 
reactions, the heating due to fission in our fit window 
(0.1 to 100 days) is solely due to spontaneous fission and 
it is close to constant during the fit window because there 
is usually one nuclide that dominates the fission heating. 
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In 85% of all cases it varies by less than a factor of two 
within the fit window, and in 99% of all cases it varies 
by less than a factor of three. Thus it is sufficient to 
report the geometric mean of the heating rate due to 
fission over the fit window. Fits to the heating rates 
over our entire parameter space are available at http: 
//stellarcollapse.org/lippunerroberts2015, 

2.6. Dominant nuclear decays 

To determine the particular nuclei that are likely to 
power kilonovae, we integrate the fractional heating con¬ 
tributions of all nuclides to find out which nuclides con¬ 
tribute most to the heating. For a single nucleosynthesis 
calculation, we know the total heating rate e(t) as a func¬ 
tion of time and we can calculate the heating rate e^(t) 
due to nuclide i as a function of time, e* (t) is calculated 
as 


Ci(t)=N A J2 (7) 


where a is an index of a reaction in the reaction net¬ 
work and it runs over the set which is the set of all 
reactions that destroy exactly one nuclide i. Na is the 
Avogadro constant in baryong -1 , A a (t) is the reaction 
rate of reaction a in s -1 , Q a is the energy released in 
reaction a in erg, and Yi(t) is the number abundance of 
nuclide i in baryon - . Note that the total heating rate 
is e(t) = where i runs over all nuclear species in 

the network. 

At any given time £, we can now calculate the fractional 
heating contribution of nuclide i as 6i(t)/e(t), which is the 
fraction of the total heating rate at time t that is solely 
due to the decay of nuclide i. These fractional heating 
contributions tell us which nuclides dominate the heating 
at a given time. To quantify which nuclides dominate the 
heating over a period of time, we define the integrated 
fractional heating contribution fi as 


fi 


1 

Inti/to 



£_dt) 

e{t) 


dint , 


(8) 


where to = 0.1 days and t\ = 100 days are the beginning 
and end of our heating rate fit window. We integrate in 
logspace to equally weigh contributions at early and late 
times. Since we know e* and e only at specific time steps 
t/c, we approximate the integral as 


h 




i 

\nti/to 


E 


^ifk) tk 


(9) 


If no tk is equal to to or ti, we add these two endpoints 
to the sum and interpolate e* and e at those points. 

Note that we calculate fi for each nuclide i in a sin¬ 
gle nucleosynthesis calculation. So we should really say 
/i(y e ,s,r), because fi will be different for the same nu¬ 
clide i in different nucleosynthesis calculations since dif¬ 
ferent amounts of nuclide i are be produced, depending 
on T e , 8, and r. To get an idea of which nuclides have 
the biggest influence on the heating rate over a range of 
T e , 8, and r, we average fi over multiple nucleosynthe¬ 
sis calculations in our parameter space. We call this the 
average integrated fractional heating contribution fi and 


calculate it as 

fi 


rni^im 


y ao) 


v e ey seSreT 


where Y, 5, and T are the sets of values of Y e , 8, and 
r, respectively, that we are averaging over, and |Y|, |<S|, 
and |T| are the cardinalities of those sets, i.e. the num¬ 
ber of elements in the sets. Note that this method of 
averaging is meaningful because we are considering the 
fractional heating contribution of nuclide i and not the 
absolute heating contribution, and furthermore, we nor¬ 
malize /i(Y e ,s,r) in the same way for each nucleosyn¬ 
thesis calculation. The final number fi that we obtain is 
a number between 0 and 1 and it tells us that nuclide i 
is responsible for this fraction of the total heating rate 
between 0.1 and 100 days averaged over a certain set of 
parameters Y e , 8, and r. Note that fi is not intended to 
be used to estimate the absolute amount of heating due 
to nuclide i, because the absolute amount of heating can 
vary greatly between the different nucleosynthesis cases 
over which we averaged to obtain fi. Rather, fi is in¬ 
tended to quantify how important different nuclides are 
in the makeup of the total radioactive heating rate over 
a wide range of possible kilonovae. This can help in¬ 
form experiments that are measuring the /3-decay prop¬ 
erties of nuclides produced in the r-process. To model 
the r-process and associated kilonovae more accurately, 
it would be more beneficial to have precise measurements 
of the /3-decay properties of nuclides that have a larger 
fi than of nuclides with smaller fi. 

Table [3] shows the 10 most dominant heating nuclides 
and their_ average integrated fractional heating contri¬ 
butions fi. The ff s are averaged over different high- 
resolution symO (symmetric fission with no free neutrons) 
runs in different Y e bins and over the entire range of en¬ 
tropies (1 fee baryon -1 < s < 100 ks baryon -1 ) and ex¬ 
pansion timescales (0.1ms < r < 500 ms). In each Y e 
bin, the nuclides are sorted with decreasing f im We only 
look at the Y e -dependence of the dominant heating nu¬ 
clides because the r-process depends very strongly on 
while it is q uite insensitive tcy entropy (e.g. Freiburghaus 


et al. 


1999 


also see Figure [l). Only the 10 mosFkiom 
mantT heating nuclides are shown here, the full table, 
and the tables of the runs with different fission reac¬ 
tions, are available at http://stellarcollapse.org/ 
lippunerroberts2015. The single most important nu¬ 
clide for heating between 0.1 and 100 days is 132 I. It 
dominates over all other nuclides by a factor of at least 3 
to 10 and it especially dominates at low initial Y e . 132 Sn 
is doubly magic (50 protons and 82 neutrons) and so it 
gets produced in high quantities in the r-process. Within 
minutes, 132 Sn decays to 132 Sb which decays to 132 Te. 
132 Te has a half-life of 3.2 days and so it decays in the 
middle of our fit window where we are looking at the 
heating contributions. But the decay of 132 Te to 132 I has 
a Q-value of only about 500 keV, while 132 I decays to 
the stable isotope 132 Xe (which is in the middle of the 
second r-process peak) with a half-life of only 2.3 hours 
and a Q-value of 3.6 MeV. Thus we get a large heating 
contribution from 132 I. 

As is to be expected, at very low Y e (between 0 and 
0.125), most of the heating comes from nuclei that form 
the second ( A ~ 130) and third (A ~ 200) r-process 
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Table 3 

Average Integrated Fractional Heating Contributions fi of the High-Resolution symO a Runs 





Ye 

Bins b 




Overall 0 

(0 < Y e < 0.5) 

0 < Y e < 0.125 

0.125 <Y e < 0.250 

0.250 < Y e < 0.375 

0.375 < Y e < 0.5 

Nuclide 

fi 

Nuclide 

fi 

Nuclide 

fi 

Nuclide 

fi 

Nuclide 

fi 

132j 

22.59% 

132j 

26.49% 

89 Sr 

9.01% 

66 Cu 

13.21% 

132j 

13.99% 

2°°Au 

4.46% 

131j 

5.52% 

72 Ga 

5.91% 

57 Ni 

10.83% 

66 Cu 

4.42% 

128 Sb 

4.26% 

128 Sb 

4.66% 

132j 

5.00% 

59 Fe 

7.47% 

89 Sr 

3.51% 

249 Bk 

4.23% 

132 Te 

3.78% 

59 Fe 

4.77% 

89 Sr 

5.21% 

57 Ni 

3.18% 

132 Te 

3.22% 

125 Sn 

3.37% 

78 As 

4.65% 

77 As 

4.79% 

59 Fe 

3.04% 

131j 

3.13% 

133j 

3.06% 

125 Sn 

3.64% 

77 Ge 

4.18% 

128 Sb 

2.67% 

252 Cf 

3.09% 

129 Sb 

2.85% 

103 Ru 

3.24% 

61 Cu 

3.20% 

131 j 

2.59% 

133j 

3.09% 

127 Sb 

2.79% 

91y 

3.08% 

62 Cu 

3.04% 

78 As 

2.27% 

202 Au 

2.89% 

140 La 

2.56% 

66 Cu 

2.97% 

56 Ni 

3.00% 

72 Ga 

2.05% 

135j 

2.65% 

129 Te 

2.25% 

112 Ag 

2.96% 

72 Ga 

2.95% 

77 Ge 

2.02% 


a Symmetric fission reactions that do not create free neutrons. 

b The fi s shown in these columns are averaged over all nucleosynthesis calculations (with different initial electron 
fractions, entropies, and expansion timescales) whose Y e falls within the Y e bin. 
c The fi s shown in this column are averaged over the entire parameter space. 


peaks. A few very heavy nuclides (A ~ 250) contribute. 
At higher Y e (between 0.125 and 0.25), the 10 most signif¬ 
icantly contributing nuclides are all in the second peak, 
since anything in the third peak and beyond is more diffi¬ 
cult to produce. The nuclides we find to be the dominant 
source of heating at low initial Y- are consistent with 
the dominant /3-decay nuclei that Metzger et al. (2010) 

k 


found. They only investigated a Y e = O.ffioutflow and 
we confirm that this result holds for a range of electron 
fractions below 0.25. 

At Y e between 0.25 and 0.375 there is a mix of sig¬ 
nificant contributes from the first (A ~ 88) and second 
peaks. There are also some iron peak elements, but most 
isotopes on the neutron-rich side of the iron peak have 
half-lives that are either too short or too long for our fit 
window. Notable exceptions are 59 Fe, 66 Ni, 67 Cu, and 
72 Ga. We do indeed see significant contributions from 
72 Ga and 59 Fe. Instead of 66 Ni, we see its /3-decay prod¬ 
uct, 66 Cu, which has a much larger Q-value (2.6 MeV 
instead of 250 keV) and a half-life of 5 minutes. 67 Cu 
does not contribute significantly because of its relatively 
low Q-value of 560 keV. Finally, at very high Y e (be¬ 
tween 0.375 and 0.5) there are significant significant con¬ 
tributes from the proton-rich side of stability around 
the iron peak. 57 Ni dominates over 56 Ni because it has 
one more neutron—thus it is a bit easier to produce in 
slightly neutron rich conditions (Y e < 0.5)—and the /3 + - 
decay Q-value of 57 Ni is a bit larger than that of 56 Ni 
(3.3MeV vs. 2.1 MeV). Both nuclides, however, have a 
half-life that is right inside our fit window, which is why 
both contribute significantly to the total heating rate. 

The cases that produce significant amounts of actinides 
also produce nuclides that undergo spontaneous fission. 
In those cases, the heating due to fission becomes domi¬ 
nant toward the end of the fit window (at about 100 days) 
but it is subdominant throughout the rest of the fit win¬ 
dow. The nuclides that contribute the most to fission 
induced heating across the entire parameter space are 
249 Bk, 252 Cf, and 241 Pu, which have average integrated 
fractional fission heating contributionsof 33%, 21%, and 
19%, resp ectively. These numbers are fi defined in Equa¬ 
tion (10) averaged over the entire parameter space, but 


the fi s of the individual nucleosynthesis calculations de¬ 
fined in Equation were calculated using only fission 
reactions in e^(t) (cl. Equation and with e(t) be¬ 
ing the total heating rate due to nssion alone. In other 
words, averaged over all runs in the entire parameter 
space and averaged in logspace over all times between 
0.1 and 100 days, 249 Bk accounts for 33% of the entire 
heating due to fission, and similarly for the other nu¬ 
clides. If /3-delayed fission were included in our reaction 
network, it would likely significantly alter the contribu¬ 
tion of fission to the heating rate at low electron fraction. 
For higher electron fractions, the neglect of beta-delayed 
fission is unlikely to be important since very little fissible 
material is produced. 

3. LIGHT CURVES 

To test how variations in the late-time nuclear heat¬ 
ing rate and composition affect possible electromagnetic 
transients associated with neutron star mergers, we cal¬ 
culate light curves using a simplified gray radiative trans¬ 
port scheme in a spherically symmetric outflow. 

3.1. Radiative transfer methods 

The ejecta is assumed to expand homologously, such 
that r = vt. The density structure of the outflow is then 
described by 


p(t,r) = p 0 (r/t) ( — 


-3 


(ii) 


SkyNet gives a heating rate e(t), which is the total 
amount of energy released per unit mass and per unit 
time due to nuclear reactions. The majority of this en¬ 
ergy is carried away by neutrinos, but some fraction, say 
/, is thermalized in the material. So fe(t) is the heating 
rate of the material due to nuclear reactions and decays. 

For homologous outflows, the velocity can be taken as 
a Lagrangian coordinate. Writing down the gray, La- 
grangian ra diative tr ansp ort e quations to first order in 
v/c (e.g. Mihalas & Weibel-Mihalas||1999), using the ve¬ 
locity as the Lagrangian coordinate, and including energy 
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release from nuclear reactions gives 
dE 
dt 
dF 
dt 
du 
dt 


v 2 t dv 

3P 


3J~ — 1 

vt 


E = —pcnF, 


— = fe + CK,(E- aT 4 ) , 
pt V ' 


( 12 ) 

(13) 

(14) 


where E is the radiation energy density, t is the time since 
merger, v is the velocity measured in units of the speed 
of light c, F is the radiation flux, p is the density given in 
Equation (fTTb , n is the opacity, a = 4cr/c is the radiation 
constant where a is the Stefan-Boltzmann constant, T is 
the temperature of the fluid, T is the Eddington factor 
(i.e. the ratio of the radiation pressure to the radiation 
energy density), u is the specific internal energy of the 
fluid, p is the fluid pressure, / is the fraction of the heat¬ 
ing rate e that is thermalized. The heating rate is not 
entirely thermalized because a large fraction of the nu¬ 
clear decay energy goes into neutrinos and gamma rays; 
neutrinos are lost from the system and gamma rays are 
only partially thermalized. To accurately calculate the 
thermalization fraction, one would need much more de¬ 
tailed information about the /3-decays than what is avail¬ 
able in REACLIB an d one would also have t o do y-ray 
transport. Following Barnes & Kasen (2013), we adopt 
/ = 0.3. 

The fluid is assumed to be a non-relativistic, non¬ 
degenerate ideal gas with molecular weight /i, so that 
the specific internal energy is u — 3T/(2/i). The gray 
transport equations are discretized in space on a stag¬ 
gered grid, with E and u defined on zone centers and 
F defined on zone edges. The resulting system of ordi¬ 
nary differential equations is then solved in time using a 
backward Euler method. Eddington factors are obtained 
by solving the static Boltzmann transport equation on 
a tangent ray grid at the beginning of a timestep. This 


method is similar to the one described in Ensman (1994), 
specialized to an homologous outflow. The zones are cho¬ 
sen to be logarithmically increasing in size moving away 
from the maximum radius. This is done to ensure that 
the radiation decoupling layer is resolved even at high 
densities. 

The density structure is assumed to be descr ibed by 
a broken power law as argued in Chevalier^ & Soker 


(1989). This choice was made mainly to facilitate com¬ 


parison with Barnes & Kasen (2013). The power law 
break and density scale are fixed to give the desired to¬ 
tal mass and total kinetic energy of the outflow. We use 

71//" - 1 O 2 1\ /T _ nnrl n i - O 1 nrVioro /i ic "1-10/3 or\nn r\ rvf 


light, for all light curve models (e.g 

(. Hotokezaka et al. 

2013a Rosswog 2013 

Foucart et al.| 

20llf— 


We note that the 
model and the one given in Equation 0 are both propor¬ 
tional to t -3 , but they have different scale factors. The 
main point of p{t) given in Equation ([l]) is to control the 
timescale over which the density changes at the time of 
nucleosynthesis (t < Is), but extrapolating this density 
to late times and assuming that it was the uniform den¬ 
sity of a ball of gas expanding with a fixed velocity would 
lead to superluminal expansion velocities in many cases. 


Equation (11) gives a much more reasonable estimate of 
the density at late times after nucleosynthesis is over. 


Calculating the exact wavelength and temperature de¬ 
pendent opacity of a mixture is extremely difficult be¬ 
cause of the large number of elements and absorption 
lines involved. Especially the lanthanide and actinide el¬ 
ement groups have very complicated line structures and 
the most sophisticated line structure and opacity cal¬ 
culations have only b een done for a few representative 
nuclides (e.g. Kasen et al. 2013). Such detailed opacity 
calculations are beyond The scope of this work and we 
use a simple prescription to compute the gray opacity k 
as a function of temperature T and composition as 

k = k f e (T) + ^ max [K Nd (T, Xi) - n Fe (T), 0], (15) 


where k p e (T) and kn d(T, Xj) are the iron and 
neodymium opacities given in Kasen et al. (2013). The 
sum runs over all lanthanide and actinide species with 
Xi being the mass fraction of a particular lanthanide or 
actinide species. We subtract the iron opacity from the 
neodymium opacity because ^Nd (T, X*) given in Kasen 
et al. (2013) is actually the opacity of a mixture contain¬ 
ing Xi neodymium and 1 — Xi iron. Our approximation 
assumes that every lanthanide and actinide contributes 
the same number of lines with the same distribution in 
energy. The opacity used in the gray calculation is taken 
to be the Planck mean opacity, which is appropriate when 
the wavelength depende nt op a city is calculated in the 
Sobolev approximation (Kasen 2015). At temperatures 
above 10 4 K, the opacities are held constant since ioniza¬ 
tion states which would have been accessed at those tem¬ 
peratures were not included in the original opacity calcu¬ 
lation and the opacities there are artificially low (Kasen 
2015). 


3.2. Dependence of kilonova light curves on the 
outflow properties 

Figure [8] shows the light curves and heating rates 
of the cases whose final abundances are shown in Fig¬ 
ure [l] In the left panel, the lanthanide-rich cases (Y e = 
0.0170.19) are about an order of magnitude dimmer than 
the lanthanide-free case (Y e = 0.25) and they peak at 
about a week instead of about a day. The effective 
temperature at peak of the lanthanide-rich cases is also 
much lower 1600 K vs. ^ 5700 K) than the tempera¬ 
ture of the lanthanide-free case. The heating rates be¬ 
tween 0.01 and 100 days, however, are almost identical 
for those three cases, so the significant differences in the 
light curves are solely due to the amount of lanthanides 
present in the ejecta and their effect on the opacity. Com¬ 
paring the cases Y e = 0.25 and Y e = 0.50, which are 
both lanthanide-free, the impact of the heating rate on 
the light curve can be seen. The heating rate is lower 
for the Y e = 0.50 case, because mostly stable nuclei are 
produced, leading to less heating. The result is that the 
light curve of the Y e = 0.50 case peaks slightly later 
(2.6days vs. 1.2 days for Y e = 0.25), is about an order of 
magnitude dimmer, and redder (spectral temperature is 
^3000K compared to ^5700K). 

In the left panel of Figure [8j the light curves for 
Y e = 0.01 and Y e = 0.19 have a small peak at very early 
times (about 0.04days). This early peak comes from 
our underestimate of the opacity at high temperatures. 
There is also a small bump at early times in the light 
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Figure 8. The light curves and heating rates of some selected nucleosynthesis calculations. Left: Y e = 0.01,0.19,0.25,0.50, s = 
10 ks baryon -1 , and r = 7.1ms. With Y e = 0.01 and Y e = 0.19 we obtain the full r-process and so the ejecta is lanthanide-rich, 
which drastically increases the opacity, resulting in a dim transient that peaks about a week after the nucleosynthesis event. This is 
in contrast to the Y e = 0.25 case, which has a very similar heating rate as the low-Ye cases, but does not produce lanthanides, and 
thus the transient is brighter and peaks earlier. The Y e = 0.50 transient is also lanthanide-free and peaks at a few days, but because 
a significant amount of stable nuclides are produced, the heating is much less, which leads to a dim transient. Right: Y e = 0.25, 
s = 1.0, 3.2,10,100 ks baryon -1 , and r = 7.1ms. As we saw in Figure fl] the s = 1.0 ks baryon -1 and s = 100 ks baryon -1 cases are 
lanthanide-rich, while s = 3.2 ks baryon -1 and s = 10 ks baryon -1 are lanthanide-free, which is clearly visible in the light curves. Even 
though s = 3.2 kB baryon -1 and s = 10 baryon -1 have essentially the same heating rate, the s = 3.2 ks baryon -1 case is significantly 
dimmer because it has a small amount of lantha nides. T he eje cta of a binary neutron star merger is expected to have entropies between 1 
and 10 ks baryon -1 (e.g. [Goriely et al.pOll |Just et al. 2015). 


curve of the Y e = 0.50 case, which is due to the behavior 
of the heating rate at early times. When determining 
the actual peak of the light curve, we neglect all peaks 
earlier than 0.5 days, unless they are more than three 
times brighter than all peaks after 0.5 days. If there are 
no peaks after 0.5 days, we pick the brightest peak that 
is more than three times brighter than the latest peak 
(which is also before 0.5days). 

The right panel of Figure [8] shows selected light curves 
with Y e = 0.25 and various initial entropies. The cases 
s = 1 kB baryon -1 and s = 100 ks baryon -1 produce 
very typical lanthanide-rich light curves, whereas s = 
10 kg baryon -1 produces a typical lanthanide-free light 
curve, and s = 3.2 ks baryon -1 produces a light curve 
that has trace amounts of lanthanides. 

In the cases where we make lanthanides at lower K e , we 
expect the peak luminosity to increase and move to ear¬ 
lier times at higher Y e when the ejecta transitions from 
lanthanide-rich to lanthanide-free, because the large con¬ 
tribution to the opacity from the lanthanides suddenly 
goes away (Kasen et al. 2013^_Tanaka & Hotokezaka 
2013). This is shown in Figure [9] When lanthanides are 
not produced, the transient generally becomes brighter, 
shorter, and bluer. We recall from Figure [3] that the 
heating rate at 1 day tends to decrease a little when lan¬ 
thanides go away. Thus the peak luminosity L p in the 
lanthanide-free cases is larger not because there is more 
heating in those cases, but because the peak occurs ear¬ 
lier (due to the smaller opacity) and the heating rate is 
always larger at earlier times than at later times. 

Looking at the time t p of the light curve in Figure 9j 
we see that the light curve peaks at about 6 days if the 
ejecta is lanthanide-rich and at about 1 day if the ejecta 
is lanthanide-free, which is consistent with earlier work 


Hotokezaka 2013). At high Y e , where we see some oscil¬ 
lations in the heating rate due to speci fic nuclides being 
produced (as explained in Section 2.4), the variation in 
the heating rate is reflected in the peak luminosity L p 
and the peak time t p . More heating results in a brighter 
transient at later times because the heating keeps the 
ejecta hotter, and thus the opacity remains high since 
more excited levels are popula ted, which i ncrea ses the 
number of optically thick lines ( Kasen et al.| 2013). Con¬ 
versely, less heating leads to a dimmer transient at earlier 
times because the ejecta is cooler and thus the opacity 
is lower. This variation is also reflected in the effective 
temperature T e # of the transient, but to a lesser degree. 
In general, lanthanide-rich transients have T e fj ~ 1600 K, 
which peaks at A ~ 1.8 pm in the infrared H and K bands. 
Lanthanide-free transients have T e ff ~ 6000 K (although 
this is a bit lower at very high Y e where the radioactive 
heating is reduced), which peaks at A ^ 480nm in the 
optical B band. 

In Figure [9j we can also clearly see that neutron- 
rich freeze-out produces very bright, very early, and 
very ultraviolet transients. The cleanest examples are 
s30t0.1 and slOOrO.l. There the luminosity ranges from 
2 x 10 41 to 10 42 ergs -1 , the effective temperature is about 
7 x 10 4 K, which peaks at A ^ 40 nm (extreme ultravio¬ 
let), and the peak occurs about an hour after the nucle- 
osynthesis even t. These results are very similar to what 


(e.g. Roberts et al. 2011 Barnes & Kasen 2013 Tanaka & 


Metzger et al.| (2015) found, however, they found peak 
effective temperatures of ^ 10 4 K, because they used 
higher opacities (k, = 30cm 2 g -1 ) since their trajecto¬ 
ries still contai ned a significa nt amount of lanthanides 
and actinides (Metzger 2015). In our case, we do not 
find significant amounts oFlanthanides or actinides if 
we obtain a strong neutron-rich freeze-out, and thus we 
get a lower opacity, which raises the effective temper- 
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Figure 9. The light curve results as a function of Y e for selected values of s and r. To show how lanthanides and neutron-rich freeze-out 
impact the lightcurve, we again show the lanthanide and actinide abundance Ma+Ac at peak and the neutron abundance X n at 10 minutes, 
which were already shown in Figure [3] Additionally, we plot the heating rate Me at peak, with M = 10 —2 Mq, the peak luminosity L p , 
peak time and the effective blackoody temperature T e ff at peak of the light curve. Note that in the neutron-rich freeze-out cases, the 
heating rate Me and the peak timescale t p go off the scales, their values are 10 44 — 10 45 ergs -1 and 15 — 30 min, respectively. As expected, 
L p follows the heating rate quite closely, except in the cases where we get a neutron-rich freeze-out. In those cases, we get a bright, very 
blue transient at early times. The exact point in Y e of the transition from a neutron-powered transient to an ordinary kilonova in this 
figure is somewhat arbitrary, since it depends on the exact method of finding the light curve peak that we choose, as explained in the text. 
Apart from the neutron-powered transients, the general trend is that we see a slightly dimmer, redder transient at later times if th e ejecta 
is lanthanide-rich, and a brighter, bluer transient at earlier times if it is lanthanide-free. This is consistent with earlier work fe.g. Barnes 
|fe Kasen|2013t . ~ 
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ature (Li & Paczynski 1998), making such a transient 
even harder to detect because it peaks deeper in the ul¬ 
traviolet. It appears that more work is needed to consis¬ 
tently model these neutron-powered transients. 

Note that the transition point in Y e in Figure [9] where 
the light curve peaks at about 1 hour to where it peaks 
at a few days is somewhat arbitrary because it depends 
on how we determine the peak in the light curve. As 
explained above, we arbitrarily decided to only consider 
peaks occurring earlier than 0.5 days if they are more 
than three times brighter than any later peaks. The jus¬ 
tification for this is that early peaks are very short and 
thus hard to detect, but in the cases where we only ob¬ 
tain a short, bright early peak, we do not want to pick 
out any later peaks that are really just the highest points 
of very shallow and long plateaus. 

We emphasize that the outflows used in this section 
were assumed to have homogeneous compositions. In 
nature, outflows from compact object mergers will have 
some spread in electron fraction and therefore have in¬ 
homogeneous compositions. Nonetheless, our simplified 
models provide guidance on the sensitivity of kilonova 
light curves to variations in the average electron frac¬ 
tion, entropy, and dynamical timescale during r-process 
nucleosynthesis. 

3.3. Mass estimates of potential kilonova observations 

Since the ejecta mass is a parameter in our simpli¬ 
fied light curve model, we can attempt to put a lower 
bound on the ejecta mass necessary to reproduce the 
possible kilonova observations. For the possible kilonova 
associated with GRB130603B, there is one observation in 
the infrared, one upper limit in the optical, and anothe r 
uppe r limit in the infrar ed at late times (Berger et a! 


2013 


Tanvir et al. 


2013). For every point in our low- 


resolution symO parameter space we compute nine light 
curves with all combinations of v/c = 0.1, 0.2, 0.3 and 
M/Mq = 0.01, 0.05,0.15. We then compute the observed 
AB magnitudes that would result from the light curve 
at the rest frame time when the observations were made, 
taking into account redshift and the actual filter response 
of the Hubble Space Telescope (HST^j Finally, we in¬ 
terpolate the resulting observed magnitudes as a func¬ 
tion of the ejecta mass to find the minimum mass that 
reproduces the observed magnitude in the near-infrared 
band (HST filter WFC3/F160W, roughly J-band in the 
rest frame) and produces an optical signal (HST filter 
WFC3/F606W, roughly B-band in rest frame) that is 
below the observed upper limits. 

We repeat the above procedure for light curves cal¬ 
culated with different heating efficiencies / (see Equa¬ 
tion (|l4|)), as the exact value of / is not known but has 
a direct influence on the brightness of the kilonova. For 
/ = 0.1, 0.3, and 0.5, we find that the minimum (over 
our entire parameter space) ejecta masses necessary to 
match the possibly observed kilonova after GRB130603B 
are 0.09, 0.03, and 0.02 solar masses, respectively. This 
is a reasonable result, as we expect the minimum mass 
necessary to produce a kilonova of equal brightness to 
decrease as the heating efficiency increases. 

If we repeat the same procedure with the poten- 

4 http://svo2.cab.inta-csic.es/svo/theory/fps3/index. 

php?mode=browse&gname=HST 


tially observed kilonova after GRB060614 (where there 
are detections in both the near-infrared (HST filter 
WFPC2/F814W, roughly R-band in rest frame) and op¬ 
tical (HST filter WFPC2/F606W, roughly V-band in rest 


(|Yang et al. 

2015 

Jin 

our light curves caJ 

cu- 


lated with / = 0.1 can match the observations, and for 
/ = 0.3 and 0.5 we require a minimum mass of 0.04 and 
0.05 solar masses, respectively. We note that a larger 
ejecta mass is needed when a larger heating efficiency is 
assumed. Because there are observations in two bands 
for GRB060614, our fits are more sensitive to the spectral 
temperature found in the light curve models than in the 
case of GRB130603B. Qualitatively, the spectral temper¬ 
ature scales inversely with the mass and prop ortionally to 
the heating efficiency (Li & Paczynski|fl998). Therefore, 
to keep a fixed spectral temperature when increasing the 
heating efficiency the total mass also must be increased. 
Our simple method for calculating the effective temper¬ 
ature is likely inadequate for detailed confrontation with 
multi-band observations, so these minimum masses are 
necessarily approximate. Another issue with this method 


does not have a homogeneous composition (e.g. Kasen 

et al. 

20151 Just et al. 2015; Metzger et al.12015; Wanajo 

et al. 

2014 

). Thererfore, to acquire more accurate esti- 
tie minimum ejected mass in these potential 

mates of tJ 


kilonova events, more sophisticated light curve model 
and hydrodynamical simu lations are required. Such an 


analysis was performed in Hotokezaka et al. (2013b) for 
GRB130603B, where they found preferred ejecta masses 
between 0.02 and 0.1 M©. 

Nevertheless, the work we have presented here will be 
very useful to estimate the masses and maybe even other 
parameters from future observations of kilonovae. With 
a sophisticated radiation transport method, one can cal¬ 
culate accurate light curves using our heating rates and 
lanthanide and actinide abundances. A consequence of 
our finding that the heating rate does not strongly de¬ 
pend on Y e in the lanthanide-rich regime (and not even 
on s and r except at very low Y e ) is that one will be able 
to quite accurately estimate the ejecta mass of future ob¬ 
served kilonovae without precisely knowing the values of 
y e , s, and r. A caveat is, however, that one needs to 
know the heating efficiency and lanthanide and actinide 
opacities well. 


4. CONCLUSIONS 

We have systematically performed nucleosynthesis cal¬ 
culations with our new nuclear reaction network SkyNet 
for a wide range of three parameters: initial electron frac¬ 
tion (0.01 < Y e < 0.5), initial entropy l/c^baryon -1 < 
s < 100 ks baryon - , and the expansion timescale 

0.1ms < r < 500 ms during nuclear burning. We ran 
the full parameter space with different fission reactions, 
but found that there were only small quantitative and 
no qualitative differences between the different fission 
reactions. We focused our attention on the amount of 
lanthanides and actinides produced and the heating rate 
between 0.1 and 100 days after the start of the nucle¬ 
osynthesis calculation, because kilonova transients are 
expected to occur in this time frame. With a spheri¬ 
cally symmetric, gray radiation transport scheme we esti- 
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mated the peak time, peak luminosity, and peak spectral 
temperature of the kilonova light curves. 

We find that the final amount of lanthanides and 
actinides depends most strongly on Y e and the ejecta 
is lanthanide-free for Y e > 0.26. However, there are 
some regions of the parameter space where the ejecta 
is lanthanide-free even for very low electron fractions. 
Specifically, at high initial entropies and small expansion 
timescales we get a neutron-rich freeze-out, which does 
not produce lanthanides, but may result in a very bright, 
very blue transient on the timescale of an hour. At small 
initial entropies and very large expansion timescales, 
there is significant late-time heating, which causes the 
composition to go back to NSE and effectively restart 
the r-process at a much higher electron fraction, which 
was raised by /3-decays. 

Since the lanthanides and actinides can increase the 
opacity of the material by a factor of ^ 100, we find 
that the peak luminosity increases by about one order 
of magnitude and the light curve peak timescale goes 
from about a week to about a day as the ejecta becomes 
lanth anide-free. This is consistent with pre v ious works 


by|Roberts et ah' (2011|); Kasen et al. 


Hotokezaka ( |20Rj ); 


( 2011 ; 
OTlpr 


Tanaka & 

rossman et al. (2014p. Tne heating 
rate at 1 day, however, remains largely unchanged and 
decreases by no more than one order of magnitude as the 
ejecta becomes lanthanide-free. Thus the increase in the 
kilonova luminosity is due to the decrease in the opacity 
when lanthanides are no longer present, which pushes 
the peak to earlier times when the heating is stronger. 
At very high Y e (> 0.4), there are large variations in 
the heating rate because single nuclides dominate the 
heating. At lower F e , the heating rate at 1 day is very 
uniform in entropy and expansion timescale because it is 
dominated by an ensemble of nuclides that average out 
to the same heating rate at 1 day even though the exact 
composition may b e ver y different. This has already been 
found in Metzger et al. (2010) and we are now confirming 
it for a larger parameter space. 

Overall, we find only weak correlation between the lan¬ 
thanide production and heating rate. Both quantities are 
quite strongly correlated with T e , but not so much with 
one another. The heating rate at 1 day is not affected 
much when the lanthanide abundance suddenly drops by 
many order of magnitude, but it slowly declines at higher 
T,. 


In Section |2.4[ we provided three linear inequalities 
involving T e , Ins, and lnr that can be used to deter¬ 
mine if the ejecta with those properties is lanthanide- 
rich or lanthanide-free. Those inequalities give the 
correct answer in 98% of all cases. We also provide 
parametric fits for the heating rates between 0.1 and 
100 days for all cases at http://stellarcollapse.org/ 
lippunerroberts2015. The mean fractional log differ¬ 
ence between the actual heating rate and our fit is no 
more than 1% in 95% of all cases. On the same website, 
we also provide an integrated fractional heating contribu¬ 
tion to give an idea of which specific nuclides contribute 
the most to the radioactive heating. 

Our nucleosynthesis code SkyNet will be released as 
free and open-source code soon. In the meantime, those 
interested can contact the authors about getting early 
access to the code. Future versions of SkyNet will also 
include neutrino interactions. Much more work needs 


to be done to accurately model the light curves of kilo- 
novae and especially to calculate the line structure and 
hence opacity of the lanthanide and actinide elements. 
We hope that our heating rate fits will be useful to other 
researchers to calculate kilonova light curves that could 
aid with detecting such events. 
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